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Summary 

Some numerical aspects of finite-difference algorithms for nonlinear multidimensional hyper- 
bolic conservation laws with stiff nonhomogenous (source) terms are discussed. If the stiffness 
is entirely dominated by the source term, a semi-implicit shock-capturing method is proposed 
provided that the Jacobian of the source terms possesses certain properties. The proposed 
semi-implicit method can be viewed as a variant of the Bussing and Murman point-implicit 
scheme with a more appropriate numerical dissipation for the computation of strong shock 
waves. However, if the stiffness is not solely dominated by the source terms, a fully implicit 
method would be a better choice. The situation is complicated by problems that are higher 
than one dimension, and the presence of stiff source terms further complicates the solution 
procedures for alternating direction implicit (ADI) methods. Several alternatives are dis- 
cussed. The primary motivation for constructing these schemes was to address thermally and 
chemically nonequilibrium flows in the hypersonic regime. Due to the unique structure of the 
eigenvalues and eigenvectors for fluid flows of this type, the computation can be simplified, thus 
providing a more efficient solution procedure than one might have anticipated. In a companion 
paper, extension and application of the semi-implicit method to three-dimensional chemically 
nonequilibrium flows in generalized coordinates will be presented. A preliminary numerical 
computation for a five-species reacting flow is shown to be oscillation-free around the shock. 


I. Motivation and Objectives 

The second-order total variation diminishing (TVD) type of schemes [1-5] originally designed 
for nonlinear scalar hyperbolic conservation laws or constant coefficient hyperbolic systems has 
been proven to be applicable to many multidimensional fluid dynamics problems (especially 
for a perfect gas). See for example, references [6-11]. In particular, the explicit and implicit 
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TVD schemes originally developed by Harten [l], Roe [12] and Davis [13], and modified and 
generalized by Yee [2] are rather simple in structure and require the least computational effort 
among the class of TVD methods. For the two-dimensional, compressible, Euler equations for 
a perfect gas, these explicit TVD methods require approximately twice the CPU time as the 
classical Lax-Wendroff and MacCormack schemes. The gain in robustness and elimination of 
spurious oscillation by TVD schemes often justifies the extra computation. On the other hand, 
straightforward extension of these methods to include stiff source terms for large systems is 
prohibitively expensive in terms of efficiency. There has been little work in efficiently general- 
izing TVD methods to large systems with stiff source terms. Here some numerical aspects and 
new procedures based on existing TVD properties, semi-implicit methods, relaxation methods, 
and the choice of the physical equation set in treating nonequilibrium flows will be discussed. 

Stiffness of Source Terms : Treating stiff problems implicitly is a common practice for the 
integration of stiff ordinary differential equations [14]. To further speed up the process, Kreiss 
[15] proposed splitting the equation into stiff and non-stiff terms and applying a special init ; al 
condition to damp out the fast wave initially. Guerra and Gustafsson [16] extended these ideas 
by applying a semi-implicit method to time-dependent hyperbolic problems with varying time 
scales. Their technique is to treat the stiff part implicitly and at the same time damp out 
the fast wave, leaving a more manageable system. Their methodology works well for nearly 
incompressible flows or transonic flows. For general hypersonic applications, the situation is 
more complicated. However, if the stiffness is dominated by the source terms and steady-state 
solutions are especially desired, one can still use this idea for hypersonic study by treating the 
source term implicitly and the rest explicitly. 

In the area of both thermally and chemically nonequilibrium flows, the terms responsible 
for the stiffness are the source terms. Kee and Miller [17] used a first-order spatial differencing 
together with a fully implicit time-marching scheme in handling combustion-related problems. 
In order to accurately capture complex flow fields, higher-order spatial differencing is more 
desirable. The most common types of higher-order spatial differencing are the second- and 
fourth-order central differences. Bussing and Murman [18J have illustrated that taking second- 
order central difference in space and treating the source term implicitly will essentially rescale 
the equation in time so that all events occur on a similar pseudo-time scale for steady-state 
chemical reaction calculations. Similarily, other researchers such as [11,19-23] have ultilized 
some kind of implicit treatment for stiff problems. 

Classical vs Modern Shock-Capturing Schemes: Most of the methods discussed above used 
classical shock-capturing methods. They only exhibit accurate results for smooth or weak 
shock solutions and are not robust enough for hypersonic computations. This is typical of clas- 
sical shock-capturing techniques that result in oscillatory solutions across the discontinuities 
and/or nonlinear instability, while the modern shock-capturing techniques such as the upwind 
and symmetric TVD schemes of Harten [l], Yee [2], Roe [3,12], Osher and Chakravarthy [4], 
and Jameson and Lax [5] do not possess these deficiencies. In particular, the symmetric TVD 
schemes appears to be more attractive than the upwind TVD schemes since the symmetric 
schemes have a smaller operation count and are less sensitive to numerical boundary-condition 
treatments [7]. An example is the current enhancement work to program LAURA by Gnoffo 
et al. [11] which has resulted in a large savings in TVD operation steps. 
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Coupling or Uncoupling the Species Equation from the Fluid Equations: In the usage of mod- 
ern shock-capturing methods like the TVD type of schemes for the chemically reacting flows, 
Carofano [10] was the first to introduce the formulism that enabled full coupling in Harten’s ex- 
plicit TVD scheme for a two-species, two-dimensional unsteady flow in Cartesian coordinates. 
To avoid solving a large system, Gnoffo and McCandless [19], uncoupled the species equations 
from the fluid dynamics equations and solved these two sets of systems of nonlinear partial dif- 
ferential equations in a time-lag fashion (loosely coupled method) by using a point-relaxation 
technique with a second-order symmetric TVD scheme of Yee [2] and an upwind TVD scheme 
of Osher and Chakravarthy [4]. Eberhardt and Brown [24] attempted to use the eigenvalues 
and eigenvectors of the fluid dynamics equations alone to obtain a “fully coupled” first-order 
explicit TVD scheme for a one-dimensional flow. The result of Eberhardt and Brown showed 
excessive smearing at the shock when compared to the true fully coupled explict TVD result. 
Their motivation for designing such a coupling procedure was to optimized the operation count 
by avoiding multiplication of large matrices. However, as will be demonstrated in the present 
work, if one makes use of the unique structure of eigenvectors and eigenvalues for fluid flow of 
this type, the fully coupled formulation can be simplified even for a large number of species, 
thus providing a more efficient solution procedure than one might have anticipated. Moreover, 
using the eigenvalues and eigenvectors for the fully coupled equation set allows one to have the 
freedom of controlling the proper amount of numerical dissipation for the individual waves [7]. 
In particular, for the two-dimensional chemically reacting flows, the number of linear waves is 
1 where ns is the number of species. Note that in order to capture contact discontinuities 
accurately, it is very important to apply the proper amount of numerical dissipation to the 
linear waves. 

Choosing Equation Set : In the present formulation for the chemically reacting flows, the global 
continuity equation is replaced by all the species continuity equations, instead of keeping 
the global continuity equation minus one species equation. Here the redunduncy in species 
equations due to elemental conservation is assumed to be eliminated in advance. Our study 
shows that the latter approach has two deficiencies, namely: (a) the solution of the absent 
species will suffer a loss of numerical accuracy due to arithmetic subtraction since it is equal to 
the total density minus the rest of the species density, and (b) the eigenvectors and its inverses, 
which are required for the proposed schemes, are slightly more complicated than in the former 
formulation. 

Proposing New Schemes : The goal of the present work is to construct second-order high- 
resolution schemes such that they are appropriate for thermally and chemically nonequilibrium 
flows in the hypersonic regime, and, at the same time, they are efficient and simple to program 
and easy to implement into existing computer codes. The approach here is to improve the 
existing technology by combining the desirable features from a few existing numerical methods 
into a single algorithm and, most importantly, take into account the specific characteristics of 
the particular physical problem. 

Two types of schemes are proposed. If the stiffness is entirely dominated by the source 
term, a semi-implicit TVD type of shock-capturing method is proposed for both time-accurate 
and steady-state calculations provided that the Jacobian of the source terms possesses certain 
properties. The proposed semi-implicit scheme can be viewed as a variant of the Bussing 
and Murman point-implicit predictor-corrector scheme [18] with a more appropriate numerical 
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dissipation in the computation of strong shock waves in the hypersonic regime and a speed up 
in convergence rate for steady-state applications. The predictor-corrector scheme of Bussing 
and Murman in turn is the expliet MacCormack scheme with the source term treated implicitly. 
The main advantage of the predictor-corrector formulation over most other one-step second- 
order TVD schemes is that, for multidimensional problems, one can incorporate the source 
terms more efficiently than in other formulisms. 

However, if the stiffness is not solely dominated by the source terms (e.g., stiffness due to 
highly irregular grid and/or viscous flows), a fully implicit method would be a better choice. 
The situation is complicated by problems that are higher than one dimension, and the pres- 
ence of stiff source terms further complicates the solution procedures for alternating direction 
implicit (ADI) methods. In fact, there seems to be no straightforward way of treating general 
stiff source terms implicitly with ADI procedures. Several alternatives will be discussed. The 
proposed fully implicit relaxation algorithm can be viewed as a variant of a fully coupled form 
of the algorithm proposed by Gnoffo and McCandless [19]. 

Due to space limitations, the material is divided into two separate papers. The mathematical 
aspects and the derviation of the numerical methods for solving the two-dimensional hypersonic 
nonequilibrium flow in cartesian coordinates are given in this paper, while the extension of the 
semi-implicit method to the three-dimensional, fully coupled chemically reacting viscous flow 
in generalized cordinates will be presented in a companion paper [25]. In the following sections, 
the proposed numerical methods for both the semi-implicit and fully implicit shock-capturing 
schemes will be breifly described. A few test examples will also be included to support the 
performance of the proposed schemes. 

II. The Governing Equations and the Basic Schemes 

In this section, the predictor-corrector form of the basic explicit TVD scheme and an example 
to illustrate the performance of the scheme are given. Also, all the necessary terms which 
are required for the basic TVD scheme for the compressible inviscid chemically reacting flow 
equations are derived. 

2.1 The Governing Equations 

Consider a two-dimensional system of nonhomogenous hyperbolic conservation laws 


dU dF(U) 3G(V) 

Tt + ~a7~ + s{u) - 

Here [/, F(U ), G(U ), and S(U) are column vectors of k components. Let A — dF/dU and 
B = dG /dU , with (a* , a£, ..., a*) and (aj, , ...,a£) being the eigenvalues of A and B. Denote 
R x and R y as the matrices whose columns are eigenvectors of A and B, and denote R~ l and 
R~ l as the inverses of R x and R y . In the case of the compressible inviscid flow equations with 
chemical reactions where the global continuity equation is replaced by all the species continuity 
equations, 
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here m — pu,n = pv , and s* represents the production of species from chemical reactions. The 
variables are the velocity components u and v, the pressure p, the energy per unit volume e, 
and the density of the tth species p x with p = P % > c% P — where ns is the number of 

species in the model and c x is the species mass fraction. Equation (2.2) assumes all of p x are 
linearly independent. If due to elemental conservation, ns should be the number of species 
after the redunduncy in the linearly dependent species is eliminated. 

The eigenvalues of A and B are 

(oi,...,o" s+3 ) = (u,...,u,u + a,u,u - a) (2.3a) 

(aj,...,a£ s+3 ) = (v,...,v,v + a,v,v - a). (2.3b) 

Here the frozen sound speed a is 

a 2 = p p + p e (H - u 2 - k 2 ) (2.4) 
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Let the grid spacing be denoted by Ax and Ay such that x — jAx and y 
the time step be denoted by At such that t — nAt where the superscript “n” 


= kAy , and let 
is used for the 
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time index and should not be confused with the n = pu in equation (2.2). Let a 1 ^,, 

R~*\ denote the quantities a l x , /?*, R " 1 evaluated at some symmetric average of and 
Uj+i'k- Similarly, let a l k ^ lJ 72*+ a, ^kl 1 denote quantities a l y , R y , R y l evaluated at some 
symmetric average of U Jf k and In the case of the chemically reacting flows, a^. + 1 , 

i? , i , i?* 1 ! , a* j, , x , and R~* x are defined similar to the ones used by Huang [26] and 
Carofano [10]. Also, the commonly used Roe’s average is no longer valid for a nonperfect gas 

[27]. 

Define 

Q j+k = R JU {Uj+i - k ~ Ui ‘ k) (210) 

as the difference (or the jump) of the characteristic variables in the locally x-direction, and 
define 

as the difference of the characteristic variables in the locally y-direction. For thermally and 
chemically nonequilibrium flows, the eigenvalues and eigenvectors have a similar structure. For 
the two-dimensional system (2.1), if nt is the number of thermal energy variables, then the 
eigenvalues in the x-direction will have (ns A nt + 1) “tz” characteristics plus u + a and u - a 
characteristics. Here the values “a” will reflect the added thermal energy variables. 

2.2 Explicit Preditor-Corrector TVD Scheme 

With the above notation, a formal extension of the scalar explicit second-order TVD method 
[2] in predictor-corrector form via the local characteristic approach for the nonlinear hyperbolic 
system (2.1) can be written as 
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Here the superscripts “(1) and (2)” designate the values of the function evaluated at the 
intermediate solutions U and The elements of the vector + ± , denoted by (</> ( + i ) s 

with / = 1 , k for a symmetric TVD scheme, are 
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(2.13a) 



( 2 - 

where a* + A are elements of (2.10) and a*. + A are eigenvalues in the i-direction. The function 
ip is 

{(J 2 + e 2 )/2e \z\<(' ^ 2 ’ 13c) 

Here t p(z) in (2.13c) is an entropy correction to \z\ where e is a small positive parameter (see 
reference [28] for a formula for e). In practical calculations, 0.05 < e < 0.2 is a commonly used 
range. Examples of the ‘limiter’ function can be expressed as 

Q l }+ , = minmod(c^_ A ,a' + A ) + minmod(a' + A ,a' + § ) - a' + A (2.13d) 

Q l j+ 1 = minmod(a*._ A ,a l j+ x,a l j+ j) (2.13e) 

Q\ + l =mi n mod[2a'_ A! 2a' + A ,2a' +§! i(a'_ A + <*'•+§)]• (2.13f) 

The minmod function of a list of arguments is equal to the smallest number in absolute value 
if the list of arguments is of the same sign, or is equal to zero if any arguments are of opposite 
sign. 

Here, the use of the nomenclature “TVD scheme” pertains to schemes satisfying the TVD 
sufficient conditions for one-dimensional nonlinear scalar hyperbolic conservation laws and/or 
constant coefficient hyperbolic systems. The second-order scheme (2.12) belongs to the class 
of symmetric TVD methods as discussed in reference [2]. By defining a more complex <j> 1 . A , 
(2.12) can be made upwind weighted and would belong to the class of upwind TVD schemes. 
The derivation is straightforward and will not be given here. 

The formulation of this scheme is purposely broken into two parts, namely, the predictor- 
corrector step of the MacCormack explicit scheme, and an appropriate conservative dissipation 
term designed in such a way that the final scheme is TVD for the one-dimensional constant 
coefficient case. This predictor-corrector TVD method is sometimes referred as the “TVD 
MacCormack” scheme. It is a slight modification of Roe’s one-step TVD Lax-Wendroff scheme. 
If one sets Q to be equation (2.13d) and xp(z) = |z|, the scalar scheme is the same as described 
in Davis [13] and Kwong [29]. The reason for choosing the predictor-corrector step instead of 
the one-step Lax-Wendroff formulation is that the predictor-corrector form provides a natural 
and efficient inclusion of the source terms for multidimensional problems. 

Another numerical aspect is that the method of extending scalar TVD schemes to nonlinear 
systems of hyperbolic conservation laws (2.1) in equation (2.12) by using the eigenvalues and 
eigenvectors of dF/dU and dG/dU , is sometimes referred to as the local characteristic ap- 
proach and is a variant of Roe’s linear Riemann solver [30], The advantages of this approach 
as opposed to Davis’s simplified approach [13] to systems are that (a) the current approach in 
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effect uses scalar schemes on each characteristic field; (b) the limiter used need not be the same 1 
for each field; e.g., one can use a more compressive limiter for the linear fields and use a less 
compressive limiter for the nonlinear fields; see references 12,7] for some numerical examples; 
(c) one can even use different schemes for different fields; (d) it is more efficient than the exact 
or approximate Riemann solvers; and (e) it provides a natural way to linearize the implicit 
TVD schemes [8,31]. For the hyperbolic system (2.1)-(2.9), the characteristic fields consist of 
two nonlinear fields “u ± a” and (ns + 1) linear fields of the same value “u.” The contact 
discontinuities are associated with the linear fields. It has been shown [2,7,9] that the two dif- 
ferent fields required different amounts of numerical dissipation (i.e., different limiters). Often 
the limiters that are designed for the linear field might give spurious oscillation or nonphysical 
solutions for the nonlinear field. 


To illustrate the performance of this scheme, figures 1 and 2 show a computation of a moving 
planar shock at an incident Mach number Ms = 2 interacting with an obstacle using the scheme 
(2.12) together with the limiter (2.13f). The calculation is for a perfect gas where S ~ 0. The 
grid-generation method used for this study is a generalized Schwartz-Christoffel transformation. 
Figure 2 shows eight sequential stages of the diffraction process. The results indicate that the 
explicit scheme (2.12) is capable of capturing even the subtle flow discontinuities such as the 
slipstreams. A detailed description of the computation is reported by Young and Yee [32]. The 
results of references [2,7,9,32] give a strong indication that algorithm (2.12) should be a good 
candidate for solving the hypersonic nonequilibrium flows, and further inspired the authors to 
improve scheme (2.12) for more efficient solution procedures. 


2.3 More Efficient Solution Procedures 


The extra computation in (2.12) compared with the classical central difference shock- 
capturing scheme such as the Lax-Wendroff method is due to the vectors (i?$) J± i. At first 
glance, the vectors a J± i and (R$)j±i involve matrix and vector multiplication of dimension 
ns + 3 for every grid point, and thus tend to discourage their adaption to problems other 
than ideal gas flows. Researchers such as Gnoffo and McCandless [19], and Eberhardt and 
Brown [24] were motivated to pursue other avenues to solve the complicated chemically react- 
ing flow problems. However, as will be demonstrated in this subsection, if one makes use of 
the unique structure of the eigenvectors and eigenvalues for fluid flow of this type, the fully 
coupled formulation can be simplified even for a large number of species, and thus be a viable 
approach. 


With straightforward manipulations, the operation for scheme (2.12) can be simplified 
tremendously. The vector a in equation (2.10), for example, can be expressed as 


8 


A p l - c l an 
A p 2 - c l aa 


A p n * r n< aa 
Uaa v bb 4 Am ) 

2 \ a a t 

- vbb 4 An 
\(aa+ u bb- ^) 

2 V a a / 


(2-14) 


with 


aa = 


E m £ 

^p,,, Ap* - p e (uAm + t>An - Ae) 


( 2 . 15 ) 

( 2 . 16 ) 


here, for example, (Ap l ) J+ i = p* +1 “Pj and if is understood that the spatial indices in (2.14)- 

(2.16) are at (jf + |,fc). Similarily, also has a very simple form. For instance, the R$ 
associated with the x-direction flux can be expressed as 
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Here the spatial indices on (2.17) are at (j + ^,Ar). As one can see, the terms in equations 
(2.14) and (2.17) due to the species equations are simple and do not require many operations. 
Therefore, the increase in the number of species equations is not as “CPU” intensive as one 
might have anticipated. 
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III. A Semi-implicit Predictor-Corrector TVD Scheme 

The explicit TVD scheme (2.12) can be used for either time-accurate or steady-state calcu- 
lation. It is second-order-accurate in time and space. However, if the source term is stiff, the 
restriction in the time step due to stability requirement is prohibitively small and (2.12) is not 
practical, especially for steady-state applications. In this section, a semi-implicit method is 
proposed. Another alternative is a fully implicit method. The basic implicit scheme and the 
related difficulty in extending the implicit method to higher dimensions with stiff source terms 
will be discussed in section IV. 


If one follows the same idea as Bussing and Murman 1 18] in treating the source term implic- 
itly, a semi-implicit predictor-corrector TVD scheme can easily be obtained. It can be written 
as a one-parameter family of time-differencing schemes for the source term; i.e., the following 
formulation includes scheme (2.12). The proposed scheme can be written as 
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Here, D is assumed to be invertible; i.e., only the type of source terms with its Jacobian 
such that D is invertible at each grid point is permissible. The parameter 8 is in the range 
0 < < 1 For 0 =£ 0, the source terms are treated implicitly. If 9 — 1, the time-differencing for 

the source term is first-order and (3.1) is suitable for steady-state calculations. But if 8 — 1/2, 
the time-differencing is second-order and (3.1) is suitable for time- accurate calculations. 

One can simplify equation (3.1) by partitioning the vectors t/, F, 0,5, D in equation (2.2) 
as follows: 
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(3.2b) 
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Here D 21 is a null matrix and D 22 is an identity matrix. With the above definitions, the 
scheme can be greatly simplified. The procedures are as follows: taking the predictor step, for 
example, one first solves for (At/ 77 )* 1 ) by 


then for AU 1 by 
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where ( r.h.s) 7 is the right-hand side of (3.3a) with all the indices “II” replaced by “I” and with 
the term At(S J j k ) n added. In other words, one only has to invert the D 11 matrix of dimension 
(ns, ns) instead of (ns + 3, ns + 3). Similarly, one can simplify the corrector step in the same 
way. Or, to explain it in another way, one solves the fluid equations 

dU JI dF n (U) dG u (U) _ 
dt + dx + dy 

explicitly, then uses the result to solve the species equations 

dU 7 dF'(U) dG*(U) = j 

dt dx dy 

explicitly, with the exception that the chemical reaction terms are treated implicitly in (3.4b). 
In the case where S 11 is not a null vector and it is not stiff, equation (3.3) is still applicable, 
except one has to add the term Af(S 7 j, ) n to the right-hand side of (3.3a). In the steady-state 
calculations where body fitted coordinates are used, one can further speed up the convergence 
rate by using a local time-stepping approach [33,34]. 

Figures 3-5 shows a preliminary test result for a three-dimensional, five-species reacting 
flow using the semi-implicit TVD scheme (3.1) together with (2.13d) as compared with an 
existing classical shock-capturing method which supplies numerical dissipation linearly [23]. 
The numerical result is shown to be oscillation-free around the shock while the time spent 
per iteration is approximately double when compared with the method used in [23]. The 


(3.4a) 

(3.4b) 
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configuration of the numerical experiment is shown in figure 3. The inflow conditions are 
p — 1 atm, the temperature T ~ 1200 K, and the Mach number M 4 for a premixed air and 
hydrogen fuel. The species considered are // 2 , 0 2 , O//, i/ 2 0 , and A 2 , with A 2 being inert. 
The so-called “global two-step chemistry model” for combustion is used. 

This result is very preliminary and the shock can be further sharpened with more com- 
pressible and/or upwind weighted limiters. The resulting research computer code for the test 
problem was developed based on an existing code developed by A. Kumar of NASA Langley, 
and K. Uenishi of Vigyan Research Inc. Their original code uses a semi-implicit classical shock- 
capturing method. A detailed description of the numerical experiments and the extension of 
(3.1) to three-dimensional, chemically reacting flows in generalized coordinates is given in the 
companion paper [25]. 


IV. A Fully Implicit TVD Method 

Another type of shock-capturing scheme that the authors are familiar with is a one-parameter 
family of explicit and implicit TVD type of schemes [2,31]. For the nonequilibrium equation 
(2.1), these schemes can be written as 


U, 


n+ 1 
3>k 


+ 1 


f _ j?n+\ 


At 


= Kk -(!-*) 


( 1 

A«A >- fc+ * 




- Ats ?,V 


k+i 


G n , , 


AtS U 


(4-la) 


Here 9 has the same meaning as before, except the scheme is now fully implicit for e^o. 

The numerical flux function k for both the upwind and symmetric TVD schemes [2,8] 

can be expressed as 




2 ( F hk + F, + ltk 




(4.1b) 


Similarly, one can define the numerical flux G j k _ h x in this manner. The elements of the 
denoted by (<f> 1 A ) 5 for a general second-order symmetric TVD scheme are 

3+2 






(4.2) 


where + I are elements of (2.10). The function \p is (2.13c). The “limiter” function Q l + ± 
can be the same as that in equations (2.13d)-(2.13f) 

The elements of the + i denoted by {4> l ^i) 1 for a second-order upwind TVD scheme, 

originally developed by Harten [l] and later modified and generalized by the first author [7,31], 
are 



U{a l ]+i ){g l }+ , + g‘) - ^(a‘ + 7y + i)«; 


j+j 


(4.3a) 


and 
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(4.3b) 


(ffj + i 0j)/ a j + i i ^ 0 


l‘ + i = 2^K+i) | Q 
Examples of the “limiter” function g l 3 can be expressed as 
9, = minrnod(a^_ ^ , a l j+ ± ) 


»>= l a i+i a i-i 


+ l a y+i a 


a i+4 +a i-* 


<* „ , 1 = ° 


J+ i 


(4.3c) 

(4.3d) 


g 3 = S * max 


0 »min(2|a , 3+ i|,5-Q ! '_i),min(|Q‘. + i| ) 2S-a^_i) ; S’ = sgn(a' + i). (4.3e) 


4.1 A Conservative Linearized Form For Steady-State Applications 

To solve for t/ n+1 in (4.1), one normally needs to solve a set of nonlinear algebraic equations 
iteratively. One way to avoid this is to linearized the implicit operator, and solve the linearized 
form by other means. Following the same procedure as in Yee [31], a conservative linearized 
form of (4.1) can be written as 


where 


/ + 


t .) + £(*;*♦» - H h-i) - 


= - — ( 

Ax V 


— TT n + l 


E=U 


u n 


Af' l+5 = i(^ + .-n'. + j)". 


The nonstandard notation 




E 


(4.4a) 

(4.4b) 

(4.4c) 

(4.4d) 

(4.4e) 


is used and 


nj +4ifc £ = /J i+i diag[^(a' + p]«7 + I i (£? i+lifc - E itk ) (4.4f) 

^ ik+ i E = diag [V’( a * ;+ 1 )] i ( E j,k+ 1 - E Jtk ). (4.4g) 

Here + l B 3i k + 1 are Jacobians of F and G evaluated at (j + 1, fc) and (j, A; + l). The value 
D j k is the Jacobian of the source term 5 evaluated at (j, A:), and Ej t k = — tfjV The 
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expression diag( 2 / ) denotes a diagonal matrix with diagonal elements z l . To compute (4.4f,g), 
a triple matrix multiplication of dimension (ns + 3, ns -f 3) has to be performed at every grid 
point. For steady-state applications, one can simplify (4.4f,g) as 


n*+± t k E - M x I{Ej+\'k - E Jik ) (4.4h) 

M y I(Ej tk ^ - E hk ) (4 4i) 

The scalar values M x and M y are 

M x = max V'(a'i) (4.4j) 

l 2 

M y = max^fl^) (4.4k) 


and 1 is the identity matrix. Note that (4.4h,i) involve scalar multiplications only. Other 
linearized treatments can be found in [31]. 

4.2 Stiff Source Terms, ADI Approaches and Relaxation Methods 

The extra stiff term D™ k on the implicit operater (4.4) complicated the solution procedures 
for the commonly used ADI procedures. Normally, if D” k are not stiff, one can reformulate 
(4.4) by an ADI procedure like the Beam and Warming [35] algorithm for an efficient solution 
process. Unfortunately, the D™ k considered here are stiff; consequently the additional higher 
order terms due to the ADI formulation can no longer be ignored. In a different context, 
Van Dalsem and Steger [36] suggested a remedy if D™ k is a diagonal matrix with identical 
diagonal elements. However, for chemically reacting flows, the matrix D™ k is full for the upper 
(ns, ns -|- 3) entities and no straightforward way of utilizing ADI approaches for nonlinear 
system cases with a general stiff source terms can be found. 

Recently, Gnoffo, McCandless and Yee [11] have successfully demonstrated the usefulness of 
a point-relaxation method on the implicit symmetric TVD scheme (4.4) for a loosely coupled 
chemically nonequilibrium flow. Here a similar point-relaxation or line-relaxation method is 
proposed for the fully coupled system (3.1). Despite the fact that a larger equation set is 
involved than in [11], the extra operations are minimized by making use of the simplification 
procedure of section II. For a point-relaxation method, the size of the matrix inversion for (4.4) 
is (ns + 3,ns + 3) as opposed to the loosely coupled method of [11], where the size is (ns, ns). 
The gain in the freedom of controlling the appropriate amount of numerical dissipation to each 
individual wave more than compensates for the extra expense. More importantly, solving the 
fully coupled system is believed to have a better convergence rate than the loosely coupled 
approach. All of the necessary terms required for the implicit scheme (4.4) are derived in 
section II. The implicit operator of (4.4a) is diagonally dominant for D J>k = 0. Therefore 
one has to make sure that the type of source term to be used will not destroy the diagonally 
dominant property which is required for some relaxation methods. 

For steady-state application, an algorithm utilizing the TVD scheme for the viscous flows 
is to difference the hyperbolic terms the same way as before, and then central-difference the 
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viscous terms. The final algorithm is the same as equation (4.4), except that the spatial central 
differencing of the viscous term is added to the right-hand side of equation (4.4). Numerical 
tests, comparison with other approaches, and recommendations will be reported in a future 
publication. 


V. Concluding Remarks 

Two numerical algorithms are proposed for hyperbolic conservation laws that are suitable 
for thermally and chemically nonequilibrium flows in the hypersonic regime. The specific prop- 
erties of the governing equations for fluid flow of this type were taken into consideration for 
more efficient solution procedures. The main areas of consideration are to minimize opera- 
tion counts, decrease truncation errors, increase the allowable time-step constraint imposed by 
the stiff source terms, and expand the shock-capturing capability beyond classical approaches. 
Details of all the considerations are described. A preliminary test problem shows certain advan- 
tages of the proposed semi-implicit high-resolution shock-capturing scheme over the classical 
ways of supplying numerical dissipation. More numerical testing and study will be pursued in 
the immediate future. 
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